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Abstract 

We demonstrate and explain that conventional finite difference schemes for 
direct numerical integration do not approximate the continuum Kardar-Parisi- 
Zhang (KPZ) equation due to microscopic roughness. The effective diffusion 
coefficient is found to be inconsistent with the nominal one. We propose a 
novel discretization in 1+1 dimensions which does not suffer from this defi- 
ciency and elucidates the reliability and limitations of direct integration ap- 
proaches. 
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The Kardar-Parisi-Zhang (KPZ) equation has been very successful in describing a class 
of dynamical self-affine interfaces [Q. Numerous simulations on discrete models for vapor 
deposition, bacterial colony growth, directed polymers, etc. show agreements with KPZ 
predictions. Being the simplest nonlinear stochastic evolution equation for interfaces, the 
KPZ equation is believed to be relevant to a large diversity of phenomena although experi- 
mental verifications has been controversial Many numerical investigations on the subject 
have concentrated on discrete models. This work focuses on another important approach, 
namely, direct numerical integration of the KPZ equation. Amar and Family first conducted 
such large-scale integrations p|. They found scaling exponents of the resulting interfaces in 
agreement with those from discrete models. This conclusion is supported subsequently by 
more accurate works indicating the validity of the KPZ approach PJ^. 

However, it has been observed that the discretized equations in the numerical integration 
of the KPZ equation admit peculiar properties not fully compatible with their continuum 
counterparts By applying Lam and Sander's inverse method 0], we will give a 

quantitative demonstration and theoretical explanation of an abnormal behavior of the dif- 
fusion coefficient. We propose a novel discretization for the numerical integration of KPZ 
interfaces in 1+1 dimensions. Our discrete equations behaves in a much more predictable 
way as proved by exact solution of its steady state properties. The results should clarify 
the reliability and limitations of conventional numerical integration techniques on the KPZ 
equation. In addition, conventional direct numerical integration schemes for the KPZ equa- 
tion are inefficient and numerically rather unstable at high nonlinearity We will give 
a quantitative evaluation of this instability. In contrast, the new discrete equations can be 
integrated substantially more efficiently with much improved stability. 

The KPZ equation gives the local rate of growth of the coarse-grained height profile 
h{x,t) of an interface at substrate position x and time t |]I|: 

^ = c + uV'h + ^{Vhr + v{x,t), (1) 

where c, p and A are the average growth rate, the diffusion coefficient and the nonlinear 
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parameter respectively. There is an implicit lower wavelength cutoff below which h is smooth. 
The noise 1] has a Gaussian distribution and mean and a correlator < ri{x,t)ri{x' ,t') >= 
2D6{x — x')6{t — t') . Most previous works on the numerical integration of the KPZ equation 
adopt the finite difference and Euler's method with the following equation 0-^: 

+ (Ao/8)(/ir+i - ht,)'] + pD.MQ, (2) 

where approximates h{iAx, nAt) with periodic boundary conditions and every is an 
independent Gaussian variable with zero mean and unit variance. The subscripted param- 
eters z/q, Aq and Dq are nominal values used in the iteration to be distinguished from the 
continuum values in the KPZ description in Eq. (|I]). Following previous works the 
spatial step size Ax is taken to be 1. Other choices will be discussed later. The implicit 
lower wavelength cutoff here is effectively 1 due to the spatial discretization. The temporal 
step size At is taken to be small enough so that decreasing its value further will not alter 
the results. 

In fact, it can be verified easily that as long as it is numerically stable, At need not 
be small before the KPZ scaling exponents can be computed accurately. This is because 
the discrete equation with finite At is by itself in the KPZ universality class similar to 
many discrete models due to symmetry considerations. The reason for taking a small At is 
to allow the discrete equation to approximate the continuum KPZ equation. However, we 
suggest that finite differencing is not a good approximation because of microscopic roughness, 
although this does not alter the scaling exponents due to universality. 

Naively assuming that h{x, t) is smooth at the lattice level, the error of discretization is 
0{Ax'^). Figure |I] shows the details of an interface generated using Eq. (Q). In all numerical 
simulations in this letter, we put = Dq = 1. Other choices can be recast into this form 
by rescaling the height and the time scales [Q. Here, we have taken Aq = 3 corresponding 
to medium nonlinearity and At = 0.01 and a smaller At gives similar results. Existence 
of microscopic roughness is evident and hence errors due to the spatial discretization are 
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uncontrolled. Therefore, there is no apriori reason why Eq. (j^) should approximate Eq. (|l]) 
unless there exist special reasons such as conservation laws as in the linear Aq = case. 

We apply the inverse method ||^ to examine interfaces generated using Eq. (|^) with the 
same parameters uo = Dq = 1, Xq = 3 and At = 0.01 on a lattice of size L = 32768. This 
approach computes all the parameters in the corresponding continuum KPZ description in 
Eq. (|l|). The method extracts the continuum parameters by requiring them, when plugged 
into Eq. (|l]), to give the best prediction on the evolution of a surface coarse-grained up to 
length / during a period r. Hence / and r are respectively the spatial and temporal resolutions 
of observation [0. Figure |^ shows the results. Continuum parameters are extracted at 
large / where finite size effects are insignificant. We obtain the unrenormalized parameter 
A ^ 3.04. The values u and D are renormalized in the same way due to a fluctuation 
dissipation theorem ||^ and we have D/z/ ~ 0.88 for all r. At large /, r controls the extent of 
renormalization since it dictates how short the wavelength of the modes should be to evolve 
fast enough to contribute to renormalization [j^. At small r corresponding to the short time 
limit in which no renormalization has taken place, we obtain the unrenormalized parameters 
D ~ 1.007 and u ~ 1.14. We thus have A ~ Aq and D c:^ Dq consistent with the validity 
of the finite difference approximation. Unfortunately, it is clear that v ^ vq. We check the 
result by calculating the ratio D jv independently from the correlation function [|l^ : 



C(r) =< \h{x + r, t) - h{x, t)f >= {D/p)r, (3) 

where the last equality is true for large r. We obtain D/z/ ^ 0.86 in reasonable agreement 
with the inverse method estimate. We have also repeated the measurements for a much 
smaller At = 0.00125. The ratio D/v estimated from both the inverse method and the 
correlation function is indistinguishable from the previous results within our statistical error 
which is less than ±0.02. Our estimates oi D /v m the range 0.86 to 0.88 are distinctly differ- 
ent from Dq/uq = 1. We conclude that the discrete equation (0) is in the KPZ universality 
class and is closely related to but does not approximate the continuum KPZ equation (^. 
This point will be explained later. 



It should be noted that decreasing Ax is not a vahd way to improve the accuracy of 
the finite difference scheme but is simply equivalent to diminishing the nonlinear parameter 
A. This is because any value of Ax can be rescaled back to 1 by the transformation x — *■ 
(Ax)~"^x, t (Ax)^^t, h {Ax)^^^'^h which leaves Eq. (Q) invariant except that A is now 
replaced by (Ax) A. In general, maintaining sufficient nonlinearity of the system is essential 
to exhibit any relevant properties of the KPZ class. It is most convenient to fix Ax = 1 and 
adjust the nonlinearity using A as in Ref. although tuning the nonlinearity with Ax 

has also been done @]. 

To further understand the anomaly, it is instructive to study the following novel dis- 
cretization which does give a correct diffusion coefficient: 



^ = i/or. + y^. + r/,(t), 
for i = 1 to L with periodic boundary conditions, where 



r? — ^i+i + hi^i — 2hi 

= (l/3)[(/ii+i - hif + {hi+i - hi){h, - hi^i) 



(4) 



(5) 



(6) 



The noise rjiit) has a Gaussian distribution and < rji{t)rjj{t') >= 2DQSijS(t — t'). Both Eqs. 
(0) and (^) could be equally valid (9(Ax^) spatial discretizations of Eq. (|l]) if the interface 
were smooth. However, neither are necessarily a good finite difference approximation of Eq. 
(HI) due to the microscopic roughness. 

The steady state properties of the new equations admit elegant exact solutions. Let 
P[h,t] be the probability distribution of the discrete interface at time t. It obeys 

the Fokker-Planck equation 



dP 
'dt 



^ d 

Y — 

^ dh- 

i=i 



-or. + ^^A p 



(7) 



which has the exact steady state solution: 



P[h] = exp 



'2Dr 



i+l 



hi 



(8) 
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The existence of this exact solution results from the vanishing of the Aq term on the RHS 
of Eq. when Eq. (|^) is applied, in complete analogy with the continuum case [|I|. The 
specific form of \l'j in Eq. (|^) is specifically chosen to allow for the cancelation and this 
property is not shared by other discretizations in general. 

We now calculate the associated continuum parameters in the KPZ description of Eq. 
(I). Assuming a large lattice, it follows from Eq. that the correlation function is 
C(r) = {Dq/uqit. Comparing with the continuum result in Eq. (|), we obtain D /u = Dq/uq. 
At the short time limit, the noise terms dominate in both Eqs. (|l|) and (H) and it is easy to see 
that the continuum short time noise parameter is D = Dq. Hence, the short time continuum 
diffusion coefficient is u = uq. To calculate A, we consider a screw boundary condition 
so that the interface has an average slope u. The steady state probability distribution 
now becomes P[h] = exp [— (z/o/2Do) Z]i(^t+i — hi — m)^]. It is then easy to show that the 
average growth velocity is v{u) = (dhi/dt) = Ao/3 + Aom^/2. The continuum nonlinear 
parameter can be calculated from A = f^(0) |T0[ and we get A = Aq. The average growth 
velocity is c = f (0) = Ao/3. Therefore, all three continuum parameters u, A and D are 
exactly the respective nominal values z/q, Aq and Dq in the new discretization. These results 
are confirmed numerically using both the inverse method and measurements of correlation 
functions. 

To gain further insights, we calculate the short time value of u directly by a novel 
analytical application of the inverse method. When calculating the continuum parameters 
c, z/ and A disregarding any higher order terms, the inverse method reduces the problem to 
the solution of a matrix equation 0. Due to an up-down symmetry of the interfaces at long 
length scales 0, the matrix is block diagonal at large spatial resolution / and the expression 
for u is simplified to 

< jdh/dtUd^h/dx^), > 

< {d^h/dx^)l > ^ ^ 

where the subscript c denotes coarse graining to length scale / before evaluation at a given 
lattice point i. The operation d/dx is usually carried out in the Fourier space. The growth 
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rate {dh/dt)c is contributed by F and \1/ in Eqs. (0)-(@) while effects of the noise vanish 
after averaging. Now the steady state distribution P[h] in Eq. (^) is invariant under the 
symmetry operation h —h, while the nonlinear term has the opposite parity from 
that of d^h/dx^. Therefore, < c{d'^h/ dx'^)c >= 0. In contrast, one can easily show that 
the remaining linear discrete diffusion term F approaches {d'^h/dx'^)c at sufficient coarse 
graining. Hence we get from Eq. (|^) that u = uq confirming our previous arguments. 

We now re-examine the conventional discretization in Eq. (|^) in light of our new results. 
In this case, the interface distribution P[h] has no simple solution in general. An exception is 
the linear Aq = case in which Eqs. @ and become identical. Then Eq. (|^) is again the 
exact steady state solution and we have u = uq, D = Dq and A = Aq = 0. A finite Aq perturbs 
the system. We found numerically at Aq = 3 and uq = Dq = 1 that the interface distribution 
P[h] is not far from that in Eq. (^. However, there is a small skewed correlation among 
the height differences of neighboring lattice points which can be exemplified by a skewness 
in the probability distribution of Fj defined in Eq. (^). As a result, the up-down symmetry 
of P[h] is broken and Eq. (g) now gives u = uq+ < ^ ci.d'^h/ dx^)c > / < {d'^h/dx^)1 >^ uq. 
The proof oi D = Dq is similar to the previous case. Our numerical results favor a slightly 
larger A than Aq instead of an equality. It can be caused by some non-trivial dependence of 
P[h] on the inclination in contrast to that for the new discretization. 

The new discretization is also a valuable tool for numerical investigations. Realizations 
of steady state interfaces could be generated directly using the exact distribution in Eq. (||). 
To simulate the dynamics, the equations can be integrated using Euler's method which has 
an error (9(At^/^) for stochastic equations 0. However, an operator splitting approach is 
far more efficient. Each iteration now consists of two half steps: 

=WL[/l",At] (10) 

=W^[/i"+l,At] (11) 

where Ul acting on {/i"}^]^ is the stochastic linear evolution operator for F and rj only in 
Eq. (^), while Un is the nonlinear evolution operator for the remaining \l/ term. By noting 
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that P[h] in Eq. (H) is also the steady state solution of both Ul and Un independently, it 
is easy to show that the continuum parameters c and A and the short time parameters u 
and D derived above remain exact. Hence the relevant dynamics is not perturbed despite 
an (9(At^/^) error in the particular realization of the interface. The linear operator Ul can 
be handled exactly and Eq. ( ]T(]| ) implies, after some algebra, 

, L , L 



hT~' = Y: Kl^h] + ^2DoAt Y: Kt,i] (12) 
i=i i=i 

where the ^"'s are independent standard Gaussian variables. The propagators and 

are computed from their Fourier coefficients: 

kl = exp(-7fcAt) (13) 

I 1 if = 

Kl = (14) 
[ {[1 - exp(-27fcAt)] /(27fcAt)}'/' if A; 

where 7^ = cos(27rfc/L) and k is an integer from to L — 1. In practice, the propagators 
are sharply peaked so that summation over only | i — j |< 5 in Eq. (O) is sufficient. The 



deterministic evolution in Eq. ([TTD can be integrated using the fourth-order Runge-Kutta 
method with an error (9(At^) which is also the overall error of our approach. 

To test the numerical stability, we simulate initially flat interfaces of size L = 128 at 
z/q = -Do = 1 with various values of Aq each for a period IOOOO/Aq. Figure ^ plots Ate 
against Aq where Ate is the critical time step just small enough to ensure stability during 
a run. For the new discrete equations integrated using the operator splitting Runge-Kutta 
approach and Euler's method, we found respectively Ate ~ Ag'''^^ and Ate ~ -^o ^ '^^ for 
Aq ^ 1. For the conventional approach of Euler's algorithm in Eq. (|]), the result is 
initially similar to that of the same method applied to the new discretization. However, 
Ate drops faster than exponentially at Aq ^ 8. We suspect that there is a critical Aq 
beyond which instability is unconditional. For Euler's method. Ate = 0.5 at Aq ^ 1 because 



of an instability of the diffusive term |11|. The above findings should be important for 



further understanding of numerical instabilities in equations for interfaces. The superior 
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performance of our discretized equations when integrated with the operator sphtting method 
is evident. The exponent —0.89 in fact represents a surprisingly good stabihty, taking into 
account that the dynamic time scale of the interface is proportional to Ag ^ for large Aq. 
Besides improved stability, it is about 30 times faster than the conventional approach for 
example at Aq = 3 for a small systematic error of less than 0.05% in the average growth 
velocity c. 

The KPZ equation can be recast using a Hope-Cole transformation into a form which 



describes directed polymers in random media. It has been reported in Refs. |12[ and ^ 
that the discretized versions of the transformed equation are numerically more stable than 
those obtained directly from the KPZ equation. Yet, our numerical work shows that their 
stability and hence the computational efficiency is only in between those of the conventional 
and our new discrete equations in 1+1 dimensions. These discretizations also suffer from the 
shortcoming that the effective diffusion coefficient is incompatible with the nominal value. 
Nevertheless, these methods work also in higher dimensions while it is not clear how our 
discretization can be generalized. 

In conclusion, we have explained that the conventional finite difference approach does 
not provide a genuine direct numerical integration of the KPZ equation since the continuum 
diffusion coefficient is incompatible with the nominal one in the discrete equations. This is 
explained by the microscopic roughness and skewness of the interface at steady state. Despite 
this anomaly, the discrete equations themselves are in the KPZ universality class and thus 
the scaling exponents measured in previous works are all valid. A novel discretization for the 
KPZ equation is studied. The continuum diffusion coefficient is hence shown analytically to 
be equal to the nominal value. However, this equality is due to subtle cancelation of terms 
in the Fokker-Planck equation and should not hold in general. This letter has focused on the 
KPZ equation in 1 + 1 dimensions. However, the conventional finite difference integration 
approach is also routinely applied to growth in higher dimensions as well as to variants of the 
KPZ equation and related problems such as the Kuramato-Sivaskinsky equation, etc. It 
should be interesting to examine limitations of those results in light of our findings. There 
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are few investigations on the generalization of the KPZ equation with higher order terms 
although they may be necessary for successful renormalization group calculations [jl3 



Our refined understanding on the integration approach should be important for generating 
such higher order terms in a controlled manner. 
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FIGURES 




FIG. 1. Snapshots of a segment of an interface generated by numerical integration. The time 
between two consecutive snapshots is 0.2 corresponding to 20 iterations. 
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FIG. 2. Inverse method results on the continuum parameters A, v and D as functions of the 
spatial and temporal resolutions I and r respectively. The dotted lines are A = 3.04 and v = 1.14. 
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FIG. 3. Largest possible time step Aic for numerical stability against Aq for Eq. (Q) integrated 
respectively by operator splitting approach (□) and Euler's method (o), and for Eq. @ integrated 
by Euler's method (x). 
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